mydata <- c(0, 10, 10, 20, 20, 20, 30, 30, 30, 30, 30, 40, 40, 50, 50, 60, 60, 70)

error <- function(pdistr, par) {
	if(par[1] <= 0 || par[2] <= 0) { return(Inf) }
	-ks.test(mydata, pdistr, par[1], par[2])$p.value
}

E <- mean(mydata)
V <- var(mydata)

resn <- optim(par=c(E,V), error, pdistr="pnorm")
resg <- optim(par=c((E^2)/V,E/V), error, pdistr="pgamma")

title <- sprintf('Error: N(%.0f,%.0f)=%.4f / G(%.0f,%.0f)=%.4f', resn$par[1], resn$par[2], resn$value, resg$par[1], resg$par[2], resg$value)
hist(mydata, freq=FALSE, main=title)
lines(dnorm(1:100, resn$par[1], resn$par[2]), col='red', lwd=2)
lines(dgamma(1:100, resg$par[1], resg$par[2]), col='blue', lwd=2)

